library(readr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

First Dataset, 4/10/24

Analyze dataset from 4/10/24, for the gels made on 4/8/24. There is an n = 3 for each thickness, with 3 gels of 200 um thickness (“thin”) and 3 gels of 1000 um thickness (“thick”). Three profiles were taken per gel, at a light intensity of PAR” 90 umol/m^2/s and a flow rate of 0.5 cm/s. Three additional profiles were taken per gel in the dark after 20 minutes of dark adaptation. The profiles were taken in the middle of the hydrogel, in the same position for all 6 profiles per gel.

Plotting profiles

files <- list.files(path = "./Polymerized040824_Measured041024", pattern = "\\.csv$", full.names = TRUE)

plot_list <- list()

for (file in files) {
  file_name <- tools::file_path_sans_ext(basename(file))
  
  file_data <- read_csv(paste(file),show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`)

  plot <- ggplot(data = file_data, aes(y = `Oxygen (µmol/L)`, x = `Depth (µm)`,
                                       color = `Condition`, shape = `Replicate`)) + 
      geom_point() + geom_line() +
      coord_flip() + scale_x_reverse() + ylim(190, 295) + 
      theme_bw() + ggtitle(file_name) +
      annotate("text", x = -1900, y = 295, label = paste("Start:", file_data$Time[1]), vjust = 1, hjust = 1, size = 3) +
      annotate("text", x = -1800, y = 295, label = paste("End:", file_data$Time[nrow(file_data)]), vjust = 1, hjust = 1, size = 3) +
    theme(legend.position = "bottom")
  
  #print(plot)
  
  plot_list[[file_name]] <- plot
}
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`
arranged_plots <- arrangeGrob(grobs = plot_list, ncol = 3)

gridExtra::grid.arrange(arranged_plots)

ggsave("Plots/Polymerized040824_Measured041024_plots.pdf", arranged_plots, width = 16, height = 9, units = "in", device = "pdf")

Notes about calculating flux

Wangpraseurt et al 2022: For each profile, the diffusive O2 flux was calculated using Fick’s first law of diffusion as described previously[31] (using a diffusion coefficient of DO2 = 2.2417 × 10−5). The effective DBL thickness δe as the intersection between the extrapolation of the linear part of the O2 slope within the boundary layer that intercepts with the O2 concentration of the bulk medium was calculated.[75]

75 - B. B. Jørgensen, N. P. Revsbech, Limnol. Oceanogr.1985, 30, 111.

31 - D. Wangpraseurt, S. You, F. Azam, G. Jacucci, O. Gaidarenko, M. Hildebrand, M. Kühl, A. G. Smith, M. P. Davey, A. Smith, D. D. Deheyn, S. Chen, S. Vignolini, Nat. Commun. 2020, 11, 1748.

“The diffusive O2 flux was calculated via Fick’s first law of diffusion for a water temperature = 25 °C and salinity = 30 using a molecular diffusion coefficient for O2 = 2.255 × 10−5 cm2 s−1(2). Gross photosynthesis was estimated via the light-dark shift method36. A flow chamber set-up provided slow laminar flow (flow rate = 0.5 cm s−1) and a fiber-optic halogen lamp (Schott KL2500, Schott, Germany) provided white light at defined levels of incident irradiance (400–700 nm) (0, 110, 220, and 1200 µmol photons m−2 s−1)2. Photosynthesis-irradiance curves were fitted to an exponential function37.”

Attempt to calculate flux

For the thick ones I am calculating the layer as at -250

# Get list of files for 1000um gels in the directory
files <- list.files(path = "./Polymerized040824_Measured041024", pattern = "^Gel040824_1000um_\\d\\.csv$")

# Create an empty dataframe to store results

results_df_1000um <- data.frame(
  File = character(),
  Scaffold_Thickness = character(),
  Replicate = character(),
  Flux = numeric(),
  END_DBL = numeric(),
  Oxygen_Surface = numeric(),
  Oxygen_End_DBL = numeric()
)

# Loop over each file
for (file in files) {
  # Read the file
  data <- read_csv(paste0("./Polymerized040824_Measured041024/",file), show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`)
  
  # Loop over each replicate
  for (replicate in unique(data$Replicate)) {
    # Filter data for the current replicate and condition
    oxygen_concentration <- data %>%
      filter(Replicate == replicate & Condition == "Light", `Depth (µm)` > -300) %>%
      pull(`Oxygen (µmol/L)`)
    depth <- data %>%
      filter(Replicate == replicate & Condition == "Light", `Depth (µm)` > -300) %>%
      pull(`Depth (µm)`)
    
    # Fit a linear regression model to estimate the slope within the diffusive boundary layer
    boundary_layer_lm <- lm(oxygen_concentration ~ depth)
    slope <- coef(boundary_layer_lm)[[2]]
    
    # Calculate flux using Fick’s first law of diffusion
    DO2 <- 2.2417e-5  # Diffusion coefficient of oxygen
    flux <- -DO2 * slope
    
    # Add the result to the dataframe
    results_df_1000um <- rbind(
      results_df_1000um,
      data.frame(
        File = file,
        Scaffold_Thickness = "1000 um",
        Replicate = replicate,
        Flux = flux,
        END_DBL = min(depth),
        Oxygen_Surface = max(oxygen_concentration),
        Oxygen_End_DBL = min(oxygen_concentration)
      )
    )
  }
}
## New names:
## New names:
## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`
results_df_1000um
##                     File Scaffold_Thickness Replicate          Flux END_DBL
## 1 Gel040824_1000um_4.csv            1000 um         1 -7.701213e-06    -250
## 2 Gel040824_1000um_4.csv            1000 um         2 -6.662409e-06    -250
## 3 Gel040824_1000um_4.csv            1000 um         3 -7.801359e-06    -250
## 4 Gel040824_1000um_5.csv            1000 um         1 -5.981432e-06    -250
## 5 Gel040824_1000um_5.csv            1000 um         2 -5.574327e-06    -250
## 6 Gel040824_1000um_5.csv            1000 um         3 -6.463321e-06    -250
## 7 Gel040824_1000um_6.csv            1000 um         1 -4.645738e-06    -250
## 8 Gel040824_1000um_6.csv            1000 um         2 -5.009059e-06    -250
## 9 Gel040824_1000um_6.csv            1000 um         3 -5.224788e-06    -250
##   Oxygen_Surface Oxygen_End_DBL
## 1        293.493        214.610
## 2        285.615        212.969
## 3        292.447        209.878
## 4        281.456        214.343
## 5        282.833        223.194
## 6        283.175        211.752
## 7        264.643        213.406
## 8        269.245        213.521
## 9        272.468        215.985

For the thinner ones, I am calculating with the boundary later at -400

# Get list of files for 200um gels in the directory
files <- list.files(path = "./Polymerized040824_Measured041024", pattern = "^Gel040824_200um_\\d\\.csv$")

# Create an empty dataframe to store results

results_df_200um <- data.frame(
  File = character(),
  Scaffold_Thickness = character(),
  Replicate = character(),
  Flux = numeric(),
  END_DBL = numeric(),
  Oxygen_Surface = numeric(),
  Oxygen_End_DBL = numeric()
)

# Loop over each file
for (file in files) {
  # Read the file
  data <- read_csv(paste0("./Polymerized040824_Measured041024/",file), show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`)
  
  # Loop over each replicate
  for (replicate in unique(data$Replicate)) {
    # Filter data for the current replicate and condition
    oxygen_concentration <- data %>%
      filter(Replicate == replicate & Condition == "Light", `Depth (µm)` > -450) %>%
      pull(`Oxygen (µmol/L)`)
    depth <- data %>%
      filter(Replicate == replicate & Condition == "Light", `Depth (µm)` > -450) %>%
      pull(`Depth (µm)`)
    
    # Fit a linear regression model to estimate the slope within the diffusive boundary layer
    boundary_layer_lm <- lm(oxygen_concentration ~ depth)
    slope <- coef(boundary_layer_lm)[[2]]
    
    # Calculate flux using Fick’s first law of diffusion
    DO2 <- 2.2417e-5  # Diffusion coefficient of oxygen
    flux <- -DO2 * slope
    
    # Add the result to the dataframe
    results_df_200um <- rbind(
      results_df_200um,
      data.frame(
        File = file,
        Scaffold_Thickness = "200 um",
        Replicate = replicate,
        Flux = flux,
        END_DBL = min(depth),
        Oxygen_Surface = max(oxygen_concentration),
        Oxygen_End_DBL = min(oxygen_concentration)
      )
    )
  }
}
## New names:
## New names:
## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`
results_df_200um
##                    File Scaffold_Thickness Replicate          Flux END_DBL
## 1 Gel040824_200um_4.csv             200 um         1 -1.848513e-06    -400
## 2 Gel040824_200um_4.csv             200 um         2 -1.999933e-06    -400
## 3 Gel040824_200um_4.csv             200 um         3 -1.833038e-06    -400
## 4 Gel040824_200um_5.csv             200 um         1 -1.499294e-06    -400
## 5 Gel040824_200um_5.csv             200 um         2 -1.423450e-06    -400
## 6 Gel040824_200um_5.csv             200 um         3 -1.639520e-06    -400
## 7 Gel040824_200um_6.csv             200 um         1 -1.447600e-06    -400
## 8 Gel040824_200um_6.csv             200 um         2 -1.388449e-06    -400
## 9 Gel040824_200um_6.csv             200 um         3 -1.516301e-06    -400
##   Oxygen_Surface Oxygen_End_DBL
## 1        235.138        204.172
## 2        237.367        204.110
## 3        235.725        204.696
## 4        230.437        204.200
## 5        229.684        204.580
## 6        232.367        204.267
## 7        231.559        206.548
## 8        229.721        205.829
## 9        230.814        206.309

Plotting

Full_results <- rbind(results_df_1000um, results_df_200um)

Full_results
##                      File Scaffold_Thickness Replicate          Flux END_DBL
## 1  Gel040824_1000um_4.csv            1000 um         1 -7.701213e-06    -250
## 2  Gel040824_1000um_4.csv            1000 um         2 -6.662409e-06    -250
## 3  Gel040824_1000um_4.csv            1000 um         3 -7.801359e-06    -250
## 4  Gel040824_1000um_5.csv            1000 um         1 -5.981432e-06    -250
## 5  Gel040824_1000um_5.csv            1000 um         2 -5.574327e-06    -250
## 6  Gel040824_1000um_5.csv            1000 um         3 -6.463321e-06    -250
## 7  Gel040824_1000um_6.csv            1000 um         1 -4.645738e-06    -250
## 8  Gel040824_1000um_6.csv            1000 um         2 -5.009059e-06    -250
## 9  Gel040824_1000um_6.csv            1000 um         3 -5.224788e-06    -250
## 10  Gel040824_200um_4.csv             200 um         1 -1.848513e-06    -400
## 11  Gel040824_200um_4.csv             200 um         2 -1.999933e-06    -400
## 12  Gel040824_200um_4.csv             200 um         3 -1.833038e-06    -400
## 13  Gel040824_200um_5.csv             200 um         1 -1.499294e-06    -400
## 14  Gel040824_200um_5.csv             200 um         2 -1.423450e-06    -400
## 15  Gel040824_200um_5.csv             200 um         3 -1.639520e-06    -400
## 16  Gel040824_200um_6.csv             200 um         1 -1.447600e-06    -400
## 17  Gel040824_200um_6.csv             200 um         2 -1.388449e-06    -400
## 18  Gel040824_200um_6.csv             200 um         3 -1.516301e-06    -400
##    Oxygen_Surface Oxygen_End_DBL
## 1         293.493        214.610
## 2         285.615        212.969
## 3         292.447        209.878
## 4         281.456        214.343
## 5         282.833        223.194
## 6         283.175        211.752
## 7         264.643        213.406
## 8         269.245        213.521
## 9         272.468        215.985
## 10        235.138        204.172
## 11        237.367        204.110
## 12        235.725        204.696
## 13        230.437        204.200
## 14        229.684        204.580
## 15        232.367        204.267
## 16        231.559        206.548
## 17        229.721        205.829
## 18        230.814        206.309
mean_flux <- Full_results %>%
  group_by(Scaffold_Thickness) %>%
  summarize(mean_flux = mean(Flux))

# Plotting the flux values for each thickness
ggplot(Full_results, aes(x = Scaffold_Thickness, y = Flux)) +
  geom_boxplot() + geom_point(aes(color = File), position = position_jitter(width = 0.05)) + 
  geom_text(data = mean_flux, aes(label = paste("Mean Flux = \n", sprintf("%.2e", mean_flux)), y = mean_flux-.0000012), vjust = 0.1, position = position_nudge(x = 0.25)) +
  labs(x = "Thickness (µm)", y = "Flux") +
  ggtitle("Flux Comparison between Thicknesses,\n same volumetric density (4/10/24)")

# Running ANOVA
anova_result <- aov(Flux ~ Scaffold_Thickness, data = Full_results)
summary(anova_result)
##                    Df    Sum Sq   Mean Sq F value   Pr(>F)    
## Scaffold_Thickness  1 9.098e-11 9.098e-11   136.6 3.01e-09 ***
## Residuals          16 1.065e-11 6.700e-13                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# What about for oxygen surface concentration?

mean_Oxygen_Surface <- Full_results %>%
  group_by(Scaffold_Thickness) %>%
  summarize(mean_Oxygen_Surface = mean(Oxygen_Surface))

# Plotting the surface O2 values for each thickness
ggplot(Full_results, aes(x = Scaffold_Thickness, y = Oxygen_Surface)) +
  geom_boxplot() + geom_point(aes(color = File), position = position_jitter(width = 0.05)) + 
  geom_text(data = mean_Oxygen_Surface, aes(label = paste("Mean Surface Oxygen = \n", round(mean_Oxygen_Surface, 1)), y = mean_Oxygen_Surface-25), vjust = 0.1, position = position_nudge(x = 0.0)) +
  labs(x = "Thickness (µm)", y = "Oxygen Surface Concentration (µmol/L)") +
  ggtitle("Oxygen Surface Concentration (µmol/L) Comparison between Thicknesses")

# Running ANOVA
anova_result <- aov(Oxygen_Surface ~ Scaffold_Thickness, data = Full_results)
summary(anova_result)
##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## Scaffold_Thickness  1  10395   10395   193.7 2.33e-10 ***
## Residuals          16    858      54                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using segmented linear regression, segmented R package

library(segmented)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
results_df_1000um_segmented <- data.frame(
  File = character(),
  Scaffold_Thickness = character(),
  Replicate = character(),
  slope = numeric(),
  Flux = numeric(),
  Breakpoint = numeric(),
  Oxygen_Surface = numeric(),
  Oxygen_breakpoint = numeric(),
  stringsAsFactors = FALSE
)

# Function to perform segmented regression and extract slope after breakpoint
perform_segmented_regression <- function(oxygen_concentration, depth) {
  # Fit segmented regression model
  fit <- segmented(lm(oxygen_concentration ~ depth), seg.Z = ~depth, psi = -250)
  
  # Extract the breakpoint(s)
  breakpoint <- fit$psi[, 2]
  
  # Find the index of the datapoint nearest to the breakpoint
  breakpoint_index <- which.min(abs(depth - breakpoint))

  # Extract the oxygen concentration at the breakpoint
  Oxygen_breakpoint <- oxygen_concentration[breakpoint_index]
  
  # Extract the oxygen concentration at depth = 0
  Oxygen_Surface <- oxygen_concentration[which(depth == 0)]
  
  # Extract the slope after the breakpoint for each segment
  lm_fit <- lm(oxygen_concentration ~ depth,
                 data = data.frame(depth = depth[breakpoint_index:length(depth)],
                                   oxygen_concentration = oxygen_concentration[breakpoint_index:length(oxygen_concentration)]))
  slope <- coef(lm_fit)["depth"]

  # Calculate fitted values for segmented regression lines
  fitted_values <- fitted(fit, newdata = data.frame(depth = depth))
  
  # Create dataframe with depth and fitted oxygen concentration values
  fitted_data <- data.frame("Depth" = depth, "Oxygen" = fitted_values)
  
  return(list(
    slope = slope,
    breakpoint = breakpoint,
    Oxygen_breakpoint = Oxygen_breakpoint,
    Oxygen_Surface = Oxygen_Surface,
    fitted_data = fitted_data
  ))
}

# Get list of CSV file names with "1000" in the title
file_list <- list.files(path = "./Polymerized040824_Measured041024", pattern = "1000.*\\.csv$", full.names = TRUE)

# Loop over files
for (file_path in file_list) {
  # Read data from file
  data <- read_csv(file_path, show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>%
    arrange(`Depth (µm)`)
  
  # Loop over replicates
  for (replicate_id in unique(data$Replicate)) {
    # Subset data for the current replicate
    replicate_data <- data %>%
      filter(Replicate == replicate_id & Condition == "Light")
    
    # Extract oxygen concentration and depth
    oxygen_concentration <- replicate_data %>% pull(`Oxygen (µmol/L)`)
    depth <- replicate_data %>% pull(`Depth (µm)`)
    
    # Perform segmented regression
    regression_results <- perform_segmented_regression(oxygen_concentration, depth)
    
    # Plot data and segmented regression lines
    p <- ggplot(replicate_data, aes(x = `Depth (µm)`, y = `Oxygen (µmol/L)`)) +
      geom_point() +
      geom_line(data = regression_results$fitted_data, aes(x = `Depth`, y = `Oxygen`), color = "red") +
      labs(title = paste("Replicate", replicate_id), x = "Depth (µm)", y = "Oxygen (µmol/L)")
    
    print(p)
    
    # Calculate flux using Fick’s first law of diffusion
    DO2 <- 2.2417e-5  # Diffusion coefficient of oxygen
    flux <- -DO2 * regression_results$slope
        
    # Store results in dataframe
    results_df_1000um_segmented <- bind_rows(results_df_1000um_segmented, data.frame(
      File = basename(file_path),
      Scaffold_Thickness = "1000 um",
      Replicate = replicate_id,
      slope = regression_results$slope,
      Flux = flux,  
      Breakpoint = regression_results$breakpoint,  
      Oxygen_Surface = regression_results$Oxygen_Surface, 
      Oxygen_breakpoint = regression_results$Oxygen_breakpoint,
      stringsAsFactors = FALSE
    ))
  }
}
## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

# Print the results dataframe
print(results_df_1000um_segmented)
##                             File Scaffold_Thickness Replicate     slope
## depth...1 Gel040824_1000um_4.csv            1000 um         1 0.4163000
## depth...2 Gel040824_1000um_4.csv            1000 um         2 0.2972034
## depth...3 Gel040824_1000um_4.csv            1000 um         3 0.3480109
## depth...4 Gel040824_1000um_5.csv            1000 um         1 0.3158300
## depth...5 Gel040824_1000um_5.csv            1000 um         2 0.2252707
## depth...6 Gel040824_1000um_5.csv            1000 um         3 0.3405020
## depth...7 Gel040824_1000um_6.csv            1000 um         1 0.3064920
## depth...8 Gel040824_1000um_6.csv            1000 um         2 0.2625240
## depth...9 Gel040824_1000um_6.csv            1000 um         3 0.2755740
##                    Flux Breakpoint Oxygen_Surface Oxygen_breakpoint
## depth...1 -9.332197e-06  -186.4981        293.493           214.610
## depth...2 -6.662409e-06  -228.7611        285.615           212.969
## depth...3 -7.801359e-06  -226.2641        292.447           209.878
## depth...4 -7.079961e-06  -207.0086        281.456           221.030
## depth...5 -5.049894e-06  -276.5534        282.833           215.546
## depth...6 -7.633033e-06  -176.8002        283.175           216.541
## depth...7 -6.870631e-06  -162.7789        264.643           219.004
## depth...8 -5.885001e-06  -178.6104        269.245           217.059
## depth...9 -6.177542e-06  -207.1288        272.468           218.221
print(results_df_1000um)
##                     File Scaffold_Thickness Replicate          Flux END_DBL
## 1 Gel040824_1000um_4.csv            1000 um         1 -7.701213e-06    -250
## 2 Gel040824_1000um_4.csv            1000 um         2 -6.662409e-06    -250
## 3 Gel040824_1000um_4.csv            1000 um         3 -7.801359e-06    -250
## 4 Gel040824_1000um_5.csv            1000 um         1 -5.981432e-06    -250
## 5 Gel040824_1000um_5.csv            1000 um         2 -5.574327e-06    -250
## 6 Gel040824_1000um_5.csv            1000 um         3 -6.463321e-06    -250
## 7 Gel040824_1000um_6.csv            1000 um         1 -4.645738e-06    -250
## 8 Gel040824_1000um_6.csv            1000 um         2 -5.009059e-06    -250
## 9 Gel040824_1000um_6.csv            1000 um         3 -5.224788e-06    -250
##   Oxygen_Surface Oxygen_End_DBL
## 1        293.493        214.610
## 2        285.615        212.969
## 3        292.447        209.878
## 4        281.456        214.343
## 5        282.833        223.194
## 6        283.175        211.752
## 7        264.643        213.406
## 8        269.245        213.521
## 9        272.468        215.985
results_df_200um_segmented <- data.frame(
  File = character(),
  Scaffold_Thickness = character(),
  Replicate = character(),
  slope = numeric(),
  Flux = numeric(),
  Breakpoint = numeric(),
  Oxygen_Surface = numeric(),
  Oxygen_breakpoint = numeric(),
  stringsAsFactors = FALSE
)

# Function to perform segmented regression and extract slope after breakpoint
perform_segmented_regression <- function(oxygen_concentration, depth) {
  # Fit segmented regression model
  fit <- segmented(lm(oxygen_concentration ~ depth), seg.Z = ~depth, psi = -400)
  
  # Extract the breakpoint(s)
  breakpoint <- fit$psi[, 2]
  
  # Find the index of the datapoint nearest to the breakpoint
  breakpoint_index <- which.min(abs(depth - breakpoint))

  # Extract the oxygen concentration at the breakpoint
  Oxygen_breakpoint <- oxygen_concentration[breakpoint_index]
  
  # Extract the oxygen concentration at depth = 0
  Oxygen_Surface <- oxygen_concentration[which(depth == 0)]
  
  # Extract the slope after the breakpoint for each segment
  lm_fit <- lm(oxygen_concentration ~ depth,
                 data = data.frame(depth = depth[breakpoint_index:length(depth)],
                                   oxygen_concentration = oxygen_concentration[breakpoint_index:length(oxygen_concentration)]))
  slope <- coef(lm_fit)["depth"]
  
  # Calculate fitted values for segmented regression lines
  fitted_values <- fitted(fit, newdata = data.frame(depth = depth))
  
  # Create dataframe with depth and fitted oxygen concentration values
  fitted_data <- data.frame("Depth" = depth, "Oxygen" = fitted_values)
  
  return(list(
    slope = slope,
    breakpoint = breakpoint,
    Oxygen_breakpoint = Oxygen_breakpoint,
    Oxygen_Surface = Oxygen_Surface,
    fitted_data = fitted_data
  ))
}

# Get list of CSV file names with "200" in the title
file_list <- list.files(path = "./Polymerized040824_Measured041024", pattern = "200.*\\.csv$", full.names = TRUE)

# Loop over files
for (file_path in file_list) {
  # Read data from file
  data <- read_csv(file_path, show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>%
    arrange(`Depth (µm)`)
  
  # Loop over replicates
  for (replicate_id in unique(data$Replicate)) {
    # Subset data for the current replicate
    replicate_data <- data %>%
      filter(Replicate == replicate_id & Condition == "Light")
    
    # Extract oxygen concentration and depth
    oxygen_concentration <- replicate_data %>% pull(`Oxygen (µmol/L)`)
    depth <- replicate_data %>% pull(`Depth (µm)`)
    
    # Perform segmented regression
    regression_results <- perform_segmented_regression(oxygen_concentration, depth)
    
    # Plot data and segmented regression lines
    p <- ggplot(replicate_data, aes(x = `Depth (µm)`, y = `Oxygen (µmol/L)`)) +
      geom_point() +
      geom_line(data = regression_results$fitted_data, aes(x = `Depth`, y = `Oxygen`), color = "red") +
      labs(title = paste("Replicate", replicate_id), x = "Depth (µm)", y = "Oxygen (µmol/L)")
    
    print(p)
    
    # Calculate flux using Fick’s first law of diffusion
    DO2 <- 2.2417e-5  # Diffusion coefficient of oxygen
    flux <- -DO2 * regression_results$slope
        
    # Store results in dataframe
    results_df_200um_segmented <- bind_rows(results_df_200um_segmented, data.frame(
      File = basename(file_path),
      Scaffold_Thickness = "200 um",
      Replicate = replicate_id,
      slope = regression_results$slope,
      Flux = flux,  
      Breakpoint = regression_results$breakpoint,  
      Oxygen_Surface = regression_results$Oxygen_Surface, 
      Oxygen_breakpoint = regression_results$Oxygen_breakpoint,
      stringsAsFactors = FALSE
    ))
  }
}
## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

# Print the results dataframe
print(results_df_200um_segmented)
##                            File Scaffold_Thickness Replicate      slope
## depth...1 Gel040824_200um_4.csv             200 um         1 0.10182714
## depth...2 Gel040824_200um_4.csv             200 um         2 0.12417314
## depth...3 Gel040824_200um_4.csv             200 um         3 0.11439429
## depth...4 Gel040824_200um_5.csv             200 um         1 0.09876400
## depth...5 Gel040824_200um_5.csv             200 um         2 0.09298171
## depth...6 Gel040824_200um_5.csv             200 um         3 0.10329143
## depth...7 Gel040824_200um_6.csv             200 um         1 0.09034171
## depth...8 Gel040824_200um_6.csv             200 um         2 0.06573738
## depth...9 Gel040824_200um_6.csv             200 um         3 0.08116286
##                    Flux Breakpoint Oxygen_Surface Oxygen_breakpoint
## depth...1 -2.282659e-06  -278.3708        235.138           206.914
## depth...2 -2.783589e-06  -269.4554        237.367           208.720
## depth...3 -2.564377e-06  -269.4545        235.725           208.342
## depth...4 -2.213993e-06  -236.2531        230.437           206.340
## depth...5 -2.084371e-06  -237.5620        229.684           207.066
## depth...6 -2.315484e-06  -266.7288        232.367           207.491
## depth...7 -2.025190e-06  -274.4906        231.559           209.808
## depth...8 -1.473635e-06  -340.0699        229.721           207.371
## depth...9 -1.819428e-06  -310.8676        230.814           207.816
print(results_df_200um)
##                    File Scaffold_Thickness Replicate          Flux END_DBL
## 1 Gel040824_200um_4.csv             200 um         1 -1.848513e-06    -400
## 2 Gel040824_200um_4.csv             200 um         2 -1.999933e-06    -400
## 3 Gel040824_200um_4.csv             200 um         3 -1.833038e-06    -400
## 4 Gel040824_200um_5.csv             200 um         1 -1.499294e-06    -400
## 5 Gel040824_200um_5.csv             200 um         2 -1.423450e-06    -400
## 6 Gel040824_200um_5.csv             200 um         3 -1.639520e-06    -400
## 7 Gel040824_200um_6.csv             200 um         1 -1.447600e-06    -400
## 8 Gel040824_200um_6.csv             200 um         2 -1.388449e-06    -400
## 9 Gel040824_200um_6.csv             200 um         3 -1.516301e-06    -400
##   Oxygen_Surface Oxygen_End_DBL
## 1        235.138        204.172
## 2        237.367        204.110
## 3        235.725        204.696
## 4        230.437        204.200
## 5        229.684        204.580
## 6        232.367        204.267
## 7        231.559        206.548
## 8        229.721        205.829
## 9        230.814        206.309
Full_results_segmented <- rbind(results_df_1000um_segmented, results_df_200um_segmented)

Full_results_segmented
##                              File Scaffold_Thickness Replicate      slope
## depth...1  Gel040824_1000um_4.csv            1000 um         1 0.41630000
## depth...2  Gel040824_1000um_4.csv            1000 um         2 0.29720343
## depth...3  Gel040824_1000um_4.csv            1000 um         3 0.34801086
## depth...4  Gel040824_1000um_5.csv            1000 um         1 0.31583000
## depth...5  Gel040824_1000um_5.csv            1000 um         2 0.22527071
## depth...6  Gel040824_1000um_5.csv            1000 um         3 0.34050200
## depth...7  Gel040824_1000um_6.csv            1000 um         1 0.30649200
## depth...8  Gel040824_1000um_6.csv            1000 um         2 0.26252400
## depth...9  Gel040824_1000um_6.csv            1000 um         3 0.27557400
## depth...11  Gel040824_200um_4.csv             200 um         1 0.10182714
## depth...21  Gel040824_200um_4.csv             200 um         2 0.12417314
## depth...31  Gel040824_200um_4.csv             200 um         3 0.11439429
## depth...41  Gel040824_200um_5.csv             200 um         1 0.09876400
## depth...51  Gel040824_200um_5.csv             200 um         2 0.09298171
## depth...61  Gel040824_200um_5.csv             200 um         3 0.10329143
## depth...71  Gel040824_200um_6.csv             200 um         1 0.09034171
## depth...81  Gel040824_200um_6.csv             200 um         2 0.06573738
## depth...91  Gel040824_200um_6.csv             200 um         3 0.08116286
##                     Flux Breakpoint Oxygen_Surface Oxygen_breakpoint
## depth...1  -9.332197e-06  -186.4981        293.493           214.610
## depth...2  -6.662409e-06  -228.7611        285.615           212.969
## depth...3  -7.801359e-06  -226.2641        292.447           209.878
## depth...4  -7.079961e-06  -207.0086        281.456           221.030
## depth...5  -5.049894e-06  -276.5534        282.833           215.546
## depth...6  -7.633033e-06  -176.8002        283.175           216.541
## depth...7  -6.870631e-06  -162.7789        264.643           219.004
## depth...8  -5.885001e-06  -178.6104        269.245           217.059
## depth...9  -6.177542e-06  -207.1288        272.468           218.221
## depth...11 -2.282659e-06  -278.3708        235.138           206.914
## depth...21 -2.783589e-06  -269.4554        237.367           208.720
## depth...31 -2.564377e-06  -269.4545        235.725           208.342
## depth...41 -2.213993e-06  -236.2531        230.437           206.340
## depth...51 -2.084371e-06  -237.5620        229.684           207.066
## depth...61 -2.315484e-06  -266.7288        232.367           207.491
## depth...71 -2.025190e-06  -274.4906        231.559           209.808
## depth...81 -1.473635e-06  -340.0699        229.721           207.371
## depth...91 -1.819428e-06  -310.8676        230.814           207.816
mean_flux <- Full_results_segmented %>%
  group_by(Scaffold_Thickness) %>%
  summarize(mean_flux = mean(Flux))

# Plotting the flux values for each thickness
flux_plot <- ggplot(Full_results_segmented, aes(x = Scaffold_Thickness, y = Flux)) +
  geom_boxplot() + geom_point(aes(color = File), position = position_jitter(width = 0.05)) + 
  geom_text(data = mean_flux, aes(label = paste("Mean Flux = \n", sprintf("%.2e", mean_flux)), y = mean_flux-.000002), vjust = 0.1, position = position_nudge(x = 0.25)) +
  labs(x = "Thickness (µm)", y = "Flux") +
  ggtitle("Flux Comparison between Thicknesses,\n same volumetric density (4/10/24)\n Segmented Regression")

flux_plot

ggsave("Plots/Polymerized040824_Measured041024_Flux_plot.pdf", flux_plot, device = "pdf")
## Saving 7 x 5 in image
# Running ANOVA
anova_result <- aov(Flux ~ Scaffold_Thickness, data = Full_results_segmented)
summary(anova_result)
##                    Df    Sum Sq   Mean Sq F value  Pr(>F)    
## Scaffold_Thickness  1 1.024e-10 1.024e-10   121.2 7.1e-09 ***
## Residuals          16 1.352e-11 8.400e-13                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# What about for oxygen surface concentration?

mean_Oxygen_Surface <- Full_results_segmented %>%
  group_by(Scaffold_Thickness) %>%
  summarize(mean_Oxygen_Surface = mean(Oxygen_Surface))

# Plotting the surface O2 values for each thickness
ggplot(Full_results_segmented, aes(x = Scaffold_Thickness, y = Oxygen_Surface)) +
  geom_boxplot() + geom_point(aes(color = File), position = position_jitter(width = 0.05)) + 
  geom_text(data = mean_Oxygen_Surface, aes(label = paste("Mean Surface Oxygen = \n", round(mean_Oxygen_Surface, 1)), y = mean_Oxygen_Surface-25), vjust = 0.1, position = position_nudge(x = 0.0)) +
  labs(x = "Thickness (µm)", y = "Oxygen Surface Concentration (µmol/L)") +
  ggtitle("Oxygen Surface Concentration (µmol/L) Comparison between Thicknesses")

# Running ANOVA
anova_result <- aov(Oxygen_Surface ~ Scaffold_Thickness, data = Full_results_segmented)
summary(anova_result)
##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## Scaffold_Thickness  1  10395   10395   193.7 2.33e-10 ***
## Residuals          16    858      54                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plotting profiles with DBL

files <- list.files(path = "./Polymerized040824_Measured041024", pattern = "\\.csv$", full.names = TRUE)

plot_list <- list()

for (file in files) {
  file_name <- tools::file_path_sans_ext(basename(file))
  
  file_data <- read_csv(paste(file),show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`)
  
  # Extract the DBL depth for the current file
  dbl_depth <- Full_results_segmented$Breakpoint[gsub("\\.csv$", "", Full_results_segmented$File) == file_name]

  plot <- ggplot(data = file_data, aes(y = `Oxygen (µmol/L)`, x = `Depth (µm)`,
                                       color = `Condition`, shape = `Replicate`)) + 
      geom_point() + geom_line() +
      coord_flip() + scale_x_reverse() + ylim(190, 295) + xlim(0, -750) + 
      theme_bw() + ggtitle(file_name) +
      annotate("text", x = -650, y = 295, label = paste("Start:", file_data$Time[1]), vjust = 1, hjust = 1, size = 3) +
      annotate("text", x = -550, y = 295, label = paste("End:", file_data$Time[nrow(file_data)]), vjust = 1, hjust = 1, size = 3) + 
      geom_vline(xintercept = dbl_depth, linetype = "dashed") +  # Add DBL depth line
      theme(legend.position = "bottom")
  
  #print(plot)
  
  plot_list[[file_name]] <- plot
}
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`
arranged_plots <- arrangeGrob(grobs = plot_list, ncol = 3)
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 120 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 120 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_line()`).
gridExtra::grid.arrange(arranged_plots)

ggsave("Plots/Polymerized040824_Measured041024_plots_DBL.pdf", arranged_plots, width = 16, height = 9, units = "in", device = "pdf")

Second Dataset, 4/12/24

Analyze dataset from 4/12/24, for the gels made on 4/10/24. There is an n = 3 for each thickness, with 3 gels of 200 um thickness (“thin”) and 3 gels of 1000 um thickness (“thick”). Three profiles were taken per gel, at a light intensity of PAR 90 umol/m^2/s and a flow rate of 0.5 cm/s. Three additional profiles were taken per gel in the dark after 20 minutes of dark adaptation. The profiles were taken in the middle of the hydrogel, in the same position for all 6 profiles per gel. For this dataset, I began taking profiles into the depth of the gel, which is shown in “positive” depths.

Plotting profiles

files <- list.files(path = "./Polymerized041024_Measured041224", pattern = "\\.csv$", full.names = TRUE)

plot_list <- list()

for (file in files) {
  file_name <- tools::file_path_sans_ext(basename(file))
  
  file_data <- read_csv(paste(file),show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>% filter(Condition == "Light" | Condition == "Dark")
  
  plot <- ggplot(data = file_data, aes(y = `Oxygen (µmol/L)`, x = `Depth (µm)`,
                                       color = `Condition`, shape = `Replicate`)) + 
      geom_point() + geom_line() +
      coord_flip() + scale_x_reverse() + xlim(1000, -1450) + ylim(140, 405) + 
      geom_vline(xintercept = 0, linetype = 3) + 
      theme_bw() + ggtitle(file_name) +
      annotate("text", x = -1400, y = 395, label = paste("Start:", file_data$Time[1]), vjust = 1, hjust = 1, size = 3) +
      annotate("text", x = -1300, y = 395, label = paste("End:", file_data$Time[nrow(file_data)]), vjust = 1, hjust = 1, size = 3) +
    theme(legend.position = "bottom")
  
  #print(plot)
  
  plot_list[[file_name]] <- plot
}
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`
arranged_plots <- arrangeGrob(grobs = plot_list, ncol = 3)
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
gridExtra::grid.arrange(arranged_plots)

ggsave("Plots/Polymerized041024_Measured041224_plots.pdf", arranged_plots, width = 16, height = 9, units = "in", device = "pdf")

Extra profiles taken not plotted above:

files <- c("./Polymerized041024_Measured041224/Gel041024_1000um_1.csv", "./Polymerized041024_Measured041224/Gel041024_200um_3.csv")

plot_list <- list()

for (file in files) {
  file_name <- tools::file_path_sans_ext(basename(file))
  
  file_data <- read_csv(paste(file),show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>% filter(Condition != "Light" & Condition != "Dark")

 plot <- ggplot(data = file_data, aes(y = `Oxygen (µmol/L)`, x = `Depth (µm)`,
                                       color = `Condition`, shape = `Replicate`)) + 
      geom_point() + geom_line() +
      coord_flip() + scale_x_reverse() + xlim(0, -1450) +  
      geom_vline(xintercept = 0, linetype = 3) + 
      theme_bw() + ggtitle(file_name) +
      annotate("text", x = -1400, y = 295, label = paste("Start:", file_data$Time[1]), vjust = 1, hjust = 1, size = 3) +
      annotate("text", x = -1300, y = 295, label = paste("End:", file_data$Time[nrow(file_data)]), vjust = 1, hjust = 1, size = 3) +
    theme(legend.position = "bottom")
  
  #print(plot)
  
  plot_list[[file_name]] <- plot
}
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`
arranged_plots <- arrangeGrob(grobs = plot_list, ncol = 3)
## Warning: Removed 30 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 30 rows containing missing values or values outside the scale range
## (`geom_line()`).
gridExtra::grid.arrange(arranged_plots)

ggsave("Plots/Polymerized041024_Measured041224_plots_extra.pdf", arranged_plots, width = 16, height = 9, units = "in", device = "pdf")

Using segmented linear regression, segmented R package

library(segmented)

# Function to perform segmented regression and extract slope after breakpoint
perform_segmented_regression <- function(oxygen_concentration, depth) {
  # Fit segmented regression model
  fit <- segmented(lm(oxygen_concentration ~ depth), seg.Z = ~depth, psi = -500)
  
  # Extract the breakpoint(s)
  breakpoint <- fit$psi[, 2]
  
  # Find the index of the datapoint nearest to the breakpoint
  breakpoint_index <- which.min(abs(depth - breakpoint))

  # Extract the oxygen concentration at the breakpoint
  Oxygen_breakpoint <- oxygen_concentration[breakpoint_index]
  
  # Extract the oxygen concentration at depth = 0
  Oxygen_Surface <- oxygen_concentration[which(depth == 0)]
  
  # Extract the slope after the breakpoint for each segment
  lm_fit <- lm(oxygen_concentration ~ depth,
                 data = data.frame(depth = depth[breakpoint_index:length(depth)],
                                   oxygen_concentration = oxygen_concentration[breakpoint_index:length(oxygen_concentration)]))
  slope <- coef(lm_fit)["depth"]
  
  # Calculate fitted values for segmented regression lines
  fitted_values <- fitted(fit, newdata = data.frame(depth = depth))
  
  # Create dataframe with depth and fitted oxygen concentration values
  fitted_data <- data.frame("Depth" = depth, "Oxygen" = fitted_values)

  return(list(
    slope = slope,
    breakpoint = breakpoint,
    Oxygen_breakpoint = Oxygen_breakpoint,
    Oxygen_Surface = Oxygen_Surface,
    fitted_data = fitted_data
  ))
}

results_df_1000um_segmented <- data.frame(
  File = character(),
  Scaffold_Thickness = character(),
  Replicate = character(),
  slope = numeric(),
  Flux = numeric(),
  Breakpoint = numeric(),
  Oxygen_Surface = numeric(),
  Oxygen_breakpoint = numeric(),
  stringsAsFactors = FALSE
)

# Get list of CSV file names with "1000" in the title
file_list <- list.files(path = "./Polymerized041024_Measured041224", pattern = "1000.*\\.csv$", full.names = TRUE)

# Loop over files
for (file_path in file_list) {
  # Read data from file
  data <- read_csv(file_path, show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>%
    arrange(`Depth (µm)`)
  
  # Loop over replicates
  for (replicate_id in unique(data$Replicate)) {
    # Subset data for the current replicate
    replicate_data <- data %>%
      filter(Replicate == replicate_id & Condition == "Light") %>%
      filter(`Depth (µm)` <= 0)
    
    # Extract oxygen concentration and depth
    oxygen_concentration <- replicate_data %>% pull(`Oxygen (µmol/L)`)
    depth <- replicate_data %>% pull(`Depth (µm)`)
    
    # Perform segmented regression
    regression_results <- perform_segmented_regression(oxygen_concentration, depth)
    
    # Plot data and segmented regression lines
    p <- ggplot(replicate_data, aes(x = `Depth (µm)`, y = `Oxygen (µmol/L)`)) +
      geom_point() +
      geom_line(data = regression_results$fitted_data, aes(x = `Depth`, y = `Oxygen`), color = "red") +
      labs(title = paste("Replicate", replicate_id), x = "Depth (µm)", y = "Oxygen (µmol/L)")
    
    print(p)
    
    # Calculate flux using Fick’s first law of diffusion
    DO2 <- 2.2417e-5  # Diffusion coefficient of oxygen
    flux <- -DO2 * regression_results$slope
        
    # Store results in dataframe
    results_df_1000um_segmented <- bind_rows(results_df_1000um_segmented, data.frame(
      File = basename(file_path),
      Scaffold_Thickness = "1000 um",
      Replicate = replicate_id,
      slope = regression_results$slope,
      Flux = flux,  
      Breakpoint = regression_results$breakpoint,  
      Oxygen_Surface = regression_results$Oxygen_Surface, 
      Oxygen_breakpoint = regression_results$Oxygen_breakpoint,
      stringsAsFactors = FALSE
    ))
  }
}
## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

# Print the results dataframe
print(results_df_1000um_segmented)
##                             File Scaffold_Thickness Replicate     slope
## depth...1 Gel041024_1000um_1.csv            1000 um         1 0.1886047
## depth...2 Gel041024_1000um_1.csv            1000 um         2 0.1865547
## depth...3 Gel041024_1000um_1.csv            1000 um         3 0.1909753
## depth...4 Gel041024_1000um_2.csv            1000 um         1 0.1564983
## depth...5 Gel041024_1000um_2.csv            1000 um         2 0.1645147
## depth...6 Gel041024_1000um_2.csv            1000 um         3 0.1813063
## depth...7 Gel041024_1000um_3.csv            1000 um         1 0.2525014
## depth...8 Gel041024_1000um_3.csv            1000 um         2 0.1921323
## depth...9 Gel041024_1000um_3.csv            1000 um         3 0.1993643
##                    Flux Breakpoint Oxygen_Surface Oxygen_breakpoint
## depth...1 -4.227951e-06  -417.4311        288.860           216.726
## depth...2 -4.181996e-06  -418.5536        289.131           216.760
## depth...3 -4.281094e-06  -389.4811        288.271           215.303
## depth...4 -3.508223e-06  -380.1215        273.166           212.490
## depth...5 -3.687925e-06  -387.7747        275.934           212.737
## depth...6 -4.064344e-06  -412.8095        280.751           211.682
## depth...7 -5.660325e-06  -321.0366        288.669           216.095
## depth...8 -4.307031e-06  -385.5178        287.524           212.723
## depth...9 -4.469150e-06  -385.2196        291.139           213.216
results_df_200um_segmented <- data.frame(
  File = character(),
  Scaffold_Thickness = character(),
  Replicate = character(),
  slope = numeric(),
  Flux = numeric(),
  Breakpoint = numeric(),
  Oxygen_Surface = numeric(),
  Oxygen_breakpoint = numeric(),
  stringsAsFactors = FALSE
)

# Get list of CSV file names with "200" in the title
file_list <- list.files(path = "./Polymerized041024_Measured041224", pattern = "200.*\\.csv$", full.names = TRUE)

# Loop over files
for (file_path in file_list) {
  # Read data from file
  data <- read_csv(file_path, show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>%
    arrange(`Depth (µm)`)
  
  # Loop over replicates
  for (replicate_id in unique(data$Replicate)) {
    # Subset data for the current replicate
    replicate_data <- data %>%
      filter(Replicate == replicate_id & Condition == "Light") %>%
      filter(`Depth (µm)` <= 0)
    
    # Extract oxygen concentration and depth
    oxygen_concentration <- replicate_data %>% pull(`Oxygen (µmol/L)`)
    depth <- replicate_data %>% pull(`Depth (µm)`)
    
    # Perform segmented regression
    regression_results <- perform_segmented_regression(oxygen_concentration, depth)
    
    # Plot data and segmented regression lines
    p <- ggplot(replicate_data, aes(x = `Depth (µm)`, y = `Oxygen (µmol/L)`)) +
      geom_point() +
      geom_line(data = regression_results$fitted_data, aes(x = `Depth`, y = `Oxygen`), color = "red") +
      labs(title = paste("Replicate", replicate_id), x = "Depth (µm)", y = "Oxygen (µmol/L)")
    
    print(p)
    
    # Calculate flux using Fick’s first law of diffusion
    DO2 <- 2.2417e-5  # Diffusion coefficient of oxygen
    flux <- -DO2 * regression_results$slope
        
    # Store results in dataframe
    results_df_200um_segmented <- bind_rows(results_df_200um_segmented, data.frame(
      File = basename(file_path),
      Scaffold_Thickness = "200 um",
      Replicate = replicate_id,
      slope = regression_results$slope,
      Flux = flux,  
      Breakpoint = regression_results$breakpoint,  
      Oxygen_Surface = regression_results$Oxygen_Surface, 
      Oxygen_breakpoint = regression_results$Oxygen_breakpoint,
      stringsAsFactors = FALSE
    ))
  }
}
## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

# Print the results dataframe
print(results_df_200um_segmented)
##                            File Scaffold_Thickness Replicate     slope
## depth...1 Gel041024_200um_1.csv             200 um         1 0.2021880
## depth...2 Gel041024_200um_1.csv             200 um         2 0.2184340
## depth...3 Gel041024_200um_1.csv             200 um         3 0.2466440
## depth...4 Gel041024_200um_2.csv             200 um         1 0.1583019
## depth...5 Gel041024_200um_2.csv             200 um         2 0.1779883
## depth...6 Gel041024_200um_2.csv             200 um         3 0.1804715
## depth...7 Gel041024_200um_3.csv             200 um         1 0.1273800
## depth...8 Gel041024_200um_3.csv             200 um         2 0.1296086
## depth...9 Gel041024_200um_3.csv             200 um         3 0.1346250
##                    Flux Breakpoint Oxygen_Surface Oxygen_breakpoint
## depth...1 -4.532448e-06  -423.0096        284.485           207.667
## depth...2 -4.896635e-06  -386.3812        289.267           205.874
## depth...3 -5.529019e-06  -385.2487        294.029           202.154
## depth...4 -3.548654e-06  -373.7823        258.055           204.006
## depth...5 -3.989964e-06  -420.7698        272.079           203.448
## depth...6 -4.045630e-06  -438.4452        279.738           201.503
## depth...7 -2.855477e-06  -276.4187        234.767           197.122
## depth...8 -2.905435e-06  -277.3252        235.705           198.117
## depth...9 -3.017889e-06  -278.9906        236.974           197.286
Full_results_segmented <- rbind(results_df_1000um_segmented, results_df_200um_segmented)

Full_results_segmented
##                              File Scaffold_Thickness Replicate     slope
## depth...1  Gel041024_1000um_1.csv            1000 um         1 0.1886047
## depth...2  Gel041024_1000um_1.csv            1000 um         2 0.1865547
## depth...3  Gel041024_1000um_1.csv            1000 um         3 0.1909753
## depth...4  Gel041024_1000um_2.csv            1000 um         1 0.1564983
## depth...5  Gel041024_1000um_2.csv            1000 um         2 0.1645147
## depth...6  Gel041024_1000um_2.csv            1000 um         3 0.1813063
## depth...7  Gel041024_1000um_3.csv            1000 um         1 0.2525014
## depth...8  Gel041024_1000um_3.csv            1000 um         2 0.1921323
## depth...9  Gel041024_1000um_3.csv            1000 um         3 0.1993643
## depth...11  Gel041024_200um_1.csv             200 um         1 0.2021880
## depth...21  Gel041024_200um_1.csv             200 um         2 0.2184340
## depth...31  Gel041024_200um_1.csv             200 um         3 0.2466440
## depth...41  Gel041024_200um_2.csv             200 um         1 0.1583019
## depth...51  Gel041024_200um_2.csv             200 um         2 0.1779883
## depth...61  Gel041024_200um_2.csv             200 um         3 0.1804715
## depth...71  Gel041024_200um_3.csv             200 um         1 0.1273800
## depth...81  Gel041024_200um_3.csv             200 um         2 0.1296086
## depth...91  Gel041024_200um_3.csv             200 um         3 0.1346250
##                     Flux Breakpoint Oxygen_Surface Oxygen_breakpoint
## depth...1  -4.227951e-06  -417.4311        288.860           216.726
## depth...2  -4.181996e-06  -418.5536        289.131           216.760
## depth...3  -4.281094e-06  -389.4811        288.271           215.303
## depth...4  -3.508223e-06  -380.1215        273.166           212.490
## depth...5  -3.687925e-06  -387.7747        275.934           212.737
## depth...6  -4.064344e-06  -412.8095        280.751           211.682
## depth...7  -5.660325e-06  -321.0366        288.669           216.095
## depth...8  -4.307031e-06  -385.5178        287.524           212.723
## depth...9  -4.469150e-06  -385.2196        291.139           213.216
## depth...11 -4.532448e-06  -423.0096        284.485           207.667
## depth...21 -4.896635e-06  -386.3812        289.267           205.874
## depth...31 -5.529019e-06  -385.2487        294.029           202.154
## depth...41 -3.548654e-06  -373.7823        258.055           204.006
## depth...51 -3.989964e-06  -420.7698        272.079           203.448
## depth...61 -4.045630e-06  -438.4452        279.738           201.503
## depth...71 -2.855477e-06  -276.4187        234.767           197.122
## depth...81 -2.905435e-06  -277.3252        235.705           198.117
## depth...91 -3.017889e-06  -278.9906        236.974           197.286
mean_flux <- Full_results_segmented %>%
  group_by(Scaffold_Thickness) %>%
  summarize(mean_flux = mean(Flux))

# Plotting the flux values for each thickness
flux_plot <- ggplot(Full_results_segmented, aes(x = Scaffold_Thickness, y = Flux)) +
  geom_boxplot() + geom_point(aes(color = File), position = position_jitter(width = 0.05)) + 
  geom_text(data = mean_flux, aes(label = paste("Mean Flux = \n", sprintf("%.2e", mean_flux)), y = mean_flux-.000002), vjust = 0.1, position = position_nudge(x = 0.25)) +
  labs(x = "Thickness (µm)", y = "Flux") +
  ggtitle("Flux Comparison between Thicknesses,\n same areal density (4/12/24)")

flux_plot

ggsave("Plots/Polymerized041024_Measured041224_Flux_plot.pdf", flux_plot, device = "pdf")
## Saving 7 x 5 in image
# Running ANOVA
anova_result <- aov(Flux ~ Scaffold_Thickness, data = Full_results_segmented)
summary(anova_result)
##                    Df    Sum Sq   Mean Sq F value Pr(>F)
## Scaffold_Thickness  1 5.230e-13 5.225e-13   0.836  0.374
## Residuals          16 9.997e-12 6.248e-13
# What about for oxygen surface concentration?

mean_Oxygen_Surface <- Full_results_segmented %>%
  group_by(Scaffold_Thickness) %>%
  summarize(mean_Oxygen_Surface = mean(Oxygen_Surface))

# Plotting the surface O2 values for each thickness
ggplot(Full_results_segmented, aes(x = Scaffold_Thickness, y = Oxygen_Surface)) +
  geom_boxplot() + geom_point(aes(color = File), position = position_jitter(width = 0.05)) + 
  geom_text(data = mean_Oxygen_Surface, aes(label = paste("Mean Surface Oxygen = \n", round(mean_Oxygen_Surface, 1)), y = mean_Oxygen_Surface-25), vjust = 0.1, position = position_nudge(x = 0.0)) +
  labs(x = "Thickness (µm)", y = "Oxygen Surface Concentration (µmol/L)") +
  ggtitle("Oxygen Surface Concentration (µmol/L) Comparison between Thicknesses")

# Running ANOVA
anova_result <- aov(Oxygen_Surface ~ Scaffold_Thickness, data = Full_results_segmented)
summary(anova_result)
##                    Df Sum Sq Mean Sq F value Pr(>F)  
## Scaffold_Thickness  1   1767    1767   5.627 0.0306 *
## Residuals          16   5025     314                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plotting profiles with DBL

files <- list.files(path = "./Polymerized041024_Measured041224", pattern = "\\.csv$", full.names = TRUE)

plot_list <- list()

for (file in files) {
  file_name <- tools::file_path_sans_ext(basename(file))
  
  file_data <- read_csv(paste(file),show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>% filter(Condition == "Light" | Condition == "Dark")
  
  # Extract the DBL depth for the current file
  dbl_depth <- Full_results_segmented$Breakpoint[gsub("\\.csv$", "", Full_results_segmented$File) == file_name]

  plot <- ggplot(data = file_data, aes(y = `Oxygen (µmol/L)`, x = `Depth (µm)`,
                                       color = `Condition`, shape = `Replicate`)) + 
      geom_point() + geom_line() +
      coord_flip() + scale_x_reverse() + xlim(0, -1000) + ylim(140, 405) +
      theme_bw() + ggtitle(file_name) +
      annotate("text", x = -900, y = 395, label = paste("Start:", file_data$Time[1]), vjust = 1, hjust = 1, size = 3) +
      annotate("text", x = -850, y = 395, label = paste("End:", file_data$Time[nrow(file_data)]), vjust = 1, hjust = 1, size = 3) +
      geom_vline(xintercept = dbl_depth, linetype = "dashed") +  # Add DBL depth line
      theme(legend.position = "bottom")
  
  #print(plot)
  
  plot_list[[file_name]] <- plot
}
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`
arranged_plots <- arrangeGrob(grobs = plot_list, ncol = 3)
## Warning: Removed 159 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 159 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 156 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 156 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 156 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 156 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 60 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 60 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 66 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 66 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 66 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 66 rows containing missing values or values outside the scale range
## (`geom_line()`).
gridExtra::grid.arrange(arranged_plots)

ggsave("Plots/Polymerized041024_Measured041224_plots_DBL.pdf", arranged_plots, width = 16, height = 9, units = "in", device = "pdf")

Third Dataset, 4/13/24

Analyze dataset from 4/13/24, for the gels (numbered 1-3) made on 4/08/24. These gels were stressed out using high light (PAR 630 umol/m^2/s) for two hours from 9AM-11AM in the morning in the incubator they were kept in and then held at ambient light (PAR = 90) until profiles were taken. There is an n = 3 for each thickness, with 3 gels of 200 um thickness (“thin”) and 3 gels of 1000 um thickness (“thick”). Three profiles were taken per gel, at a light intensity of PAR 90 umol/m^2/s and a flow rate of 0.5 cm/s. Three additional profiles were taken per gel in high light (PAR 630 umol/m^2/s) after 5-10 minutes of adaptation to equilibrium based on the oxygen surface concentration no longer increasing and reaching steady state. The profiles were taken in the middle of the hydrogel, in the same position for all 6 profiles per gel. Profiles were taken partially into the depth of the gel, which is shown in “positive” depths.

Plotting profiles

files <- list.files(path = "./Polymerized040824_Measured041324", pattern = "\\.csv$", full.names = TRUE)

plot_list <- list()

for (file in files) {
  file_name <- tools::file_path_sans_ext(basename(file))
  
  file_data <- read_csv(paste(file),show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) #%>% filter(Condition == "Light" | Condition == "Dark")
  
  plot <- ggplot(data = file_data, aes(y = `Oxygen (µmol/L)`, x = `Depth (µm)`,
                                       color = `Condition`, shape = `Replicate`)) + 
      geom_point() + geom_line() +
      coord_flip() + scale_x_reverse() + xlim(1000, -1450) + ylim(200, 900) + 
      geom_vline(xintercept = 0, linetype = 3) + 
      theme_bw() + ggtitle(file_name) +
      annotate("text", x = -1400, y = 395, label = paste("Start:", file_data$Time[1]), vjust = 1, hjust = 1, size = 3) +
      annotate("text", x = -1300, y = 395, label = paste("End:", file_data$Time[nrow(file_data)]), vjust = 1, hjust = 1, size = 3) +
    theme(legend.position = "bottom")
  
  #print(plot)
  
  plot_list[[file_name]] <- plot
}
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`
arranged_plots <- arrangeGrob(grobs = plot_list, ncol = 3)

gridExtra::grid.arrange(arranged_plots)

ggsave("Plots/Polymerized040824_Measured041324_plots.pdf", arranged_plots, width = 16, height = 9, units = "in", device = "pdf")

Using segmented linear regression, segmented R package

library(segmented)

# Function to perform segmented regression and extract slope after breakpoint
perform_segmented_regression <- function(oxygen_concentration, depth) {
  # Fit segmented regression model
  fit <- segmented(lm(oxygen_concentration ~ depth), seg.Z = ~depth, psi = -250)
  
  # Extract the breakpoint(s)
  breakpoint <- fit$psi[, 2]
  
  # Find the index of the datapoint nearest to the breakpoint
  breakpoint_index <- which.min(abs(depth - breakpoint))

  # Extract the oxygen concentration at the breakpoint
  Oxygen_breakpoint <- oxygen_concentration[breakpoint_index]
  
  # Extract the oxygen concentration at depth = 0
  Oxygen_Surface <- oxygen_concentration[which(depth == 0)]
  
  # Extract the slope after the breakpoint for each segment
  lm_fit <- lm(oxygen_concentration ~ depth,
                 data = data.frame(depth = depth[breakpoint_index:length(depth)],
                                   oxygen_concentration = oxygen_concentration[breakpoint_index:length(oxygen_concentration)]))
  slope <- coef(lm_fit)["depth"]
  
  # Calculate fitted values for segmented regression lines
  fitted_values <- fitted(fit, newdata = data.frame(depth = depth))
  
  # Create dataframe with depth and fitted oxygen concentration values
  fitted_data <- data.frame("Depth" = depth, "Oxygen" = fitted_values)

  return(list(
    slope = slope,
    breakpoint = breakpoint,
    Oxygen_breakpoint = Oxygen_breakpoint,
    Oxygen_Surface = Oxygen_Surface,
    fitted_data = fitted_data
  ))
}

results_df_1000um_segmented <- data.frame(
  File = character(),
  Scaffold_Thickness = character(),
  Replicate = character(),
  slope = numeric(),
  Flux = numeric(),
  Breakpoint = numeric(),
  Oxygen_Surface = numeric(),
  Oxygen_breakpoint = numeric(),
  stringsAsFactors = FALSE
)

# Get list of CSV file names with "1000" in the title
file_list <- list.files(path = "./Polymerized040824_Measured041324", pattern = "1000.*\\.csv$", full.names = TRUE)

# Loop over files
for (file_path in file_list) {
  # Read data from file
  data <- read_csv(file_path, show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>%
    arrange(`Depth (µm)`)
  
  # Loop over replicates
  for (replicate_id in unique(data$Replicate)) {
    # Subset data for the current replicate
    replicate_data <- data %>%
      filter(Replicate == replicate_id & Condition == "Light")  %>%
      filter(`Depth (µm)` <= 0)
    
    # Extract oxygen concentration and depth
    oxygen_concentration <- replicate_data %>% pull(`Oxygen (µmol/L)`)
    depth <- replicate_data %>% pull(`Depth (µm)`)
    
    # Perform segmented regression
    regression_results <- perform_segmented_regression(oxygen_concentration, depth)
    
    # Plot data and segmented regression lines
    p <- ggplot(replicate_data, aes(x = `Depth (µm)`, y = `Oxygen (µmol/L)`)) +
      geom_point() +
      geom_line(data = regression_results$fitted_data, aes(x = `Depth`, y = `Oxygen`), color = "red") +
      labs(title = paste("Replicate", replicate_id), x = "Depth (µm)", y = "Oxygen (µmol/L)")
    
    print(p)
    
    # Calculate flux using Fick’s first law of diffusion
    DO2 <- 2.2417e-5  # Diffusion coefficient of oxygen
    flux <- -DO2 * regression_results$slope
        
    # Store results in dataframe
    results_df_1000um_segmented <- bind_rows(results_df_1000um_segmented, data.frame(
      File = basename(file_path),
      Scaffold_Thickness = "1000 um",
      Replicate = replicate_id,
      slope = regression_results$slope,
      Flux = flux,  
      Breakpoint = regression_results$breakpoint,  
      Oxygen_Surface = regression_results$Oxygen_Surface, 
      Oxygen_breakpoint = regression_results$Oxygen_breakpoint,
      stringsAsFactors = FALSE
    ))
  }
}
## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

# Print the results dataframe
print(results_df_1000um_segmented)
##                             File Scaffold_Thickness Replicate     slope
## depth...1 Gel040824_1000um_1.csv            1000 um         1 0.5821080
## depth...2 Gel040824_1000um_1.csv            1000 um         2 0.5695080
## depth...3 Gel040824_1000um_1.csv            1000 um         3 0.5511640
## depth...4 Gel040824_1000um_2.csv            1000 um         1 0.5025440
## depth...5 Gel040824_1000um_2.csv            1000 um         2 0.4722500
## depth...6 Gel040824_1000um_2.csv            1000 um         3 0.4822000
## depth...7 Gel040824_1000um_3.csv            1000 um         1 0.3841909
## depth...8 Gel040824_1000um_3.csv            1000 um         2 0.4481669
## depth...9 Gel040824_1000um_3.csv            1000 um         3 0.4936017
##                    Flux Breakpoint Oxygen_Surface Oxygen_breakpoint
## depth...1 -1.304912e-05  -222.7569        350.297           235.901
## depth...2 -1.276666e-05  -214.8453        342.694           230.712
## depth...3 -1.235544e-05  -239.9065        360.905           227.053
## depth...4 -1.126553e-05  -223.8684        327.998           231.279
## depth...5 -1.058643e-05  -217.2484        317.507           225.020
## depth...6 -1.080948e-05  -214.3318        318.232           223.906
## depth...7 -8.612406e-06  -234.9614        312.305           218.707
## depth...8 -1.004656e-05  -231.9015        328.562           218.860
## depth...9 -1.106507e-05  -229.3515        341.538           220.488
results_df_200um_segmented <- data.frame(
  File = character(),
  Scaffold_Thickness = character(),
  Replicate = character(),
  slope = numeric(),
  Flux = numeric(),
  Breakpoint = numeric(),
  Oxygen_Surface = numeric(),
  Oxygen_breakpoint = numeric(),
  stringsAsFactors = FALSE
)

# Get list of CSV file names with "200" in the title
file_list <- list.files(path = "./Polymerized040824_Measured041324", pattern = "200.*\\.csv$", full.names = TRUE)

# Loop over files
for (file_path in file_list) {
  # Read data from file
  data <- read_csv(file_path, show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>%
    arrange(`Depth (µm)`)
  
  # Loop over replicates
  for (replicate_id in unique(data$Replicate)) {
    # Subset data for the current replicate
    replicate_data <- data %>%
      filter(Replicate == replicate_id & Condition == "Light") %>%
      filter(`Depth (µm)` <= 0)
    
    # Extract oxygen concentration and depth
    oxygen_concentration <- replicate_data %>% pull(`Oxygen (µmol/L)`)
    depth <- replicate_data %>% pull(`Depth (µm)`)
    
    # Perform segmented regression
    regression_results <- perform_segmented_regression(oxygen_concentration, depth)
    
    # Plot data and segmented regression lines
    p <- ggplot(replicate_data, aes(x = `Depth (µm)`, y = `Oxygen (µmol/L)`)) +
      geom_point() +
      geom_line(data = regression_results$fitted_data, aes(x = `Depth`, y = `Oxygen`), color = "red") +
      labs(title = paste("Replicate", replicate_id), x = "Depth (µm)", y = "Oxygen (µmol/L)")
    
    print(p)
    
    # Calculate flux using Fick’s first law of diffusion
    DO2 <- 2.2417e-5  # Diffusion coefficient of oxygen
    flux <- -DO2 * regression_results$slope
        
    # Store results in dataframe
    results_df_200um_segmented <- bind_rows(results_df_200um_segmented, data.frame(
      File = basename(file_path),
      Scaffold_Thickness = "200 um",
      Replicate = replicate_id,
      slope = regression_results$slope,
      Flux = flux,  
      Breakpoint = regression_results$breakpoint,  
      Oxygen_Surface = regression_results$Oxygen_Surface, 
      Oxygen_breakpoint = regression_results$Oxygen_breakpoint,
      stringsAsFactors = FALSE
    ))
  }
}
## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

## New names:
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`

# Print the results dataframe
print(results_df_200um_segmented)
##                            File Scaffold_Thickness Replicate     slope
## depth...1 Gel040824_200um_1.csv             200 um         1 0.2072080
## depth...2 Gel040824_200um_1.csv             200 um         2 0.2168434
## depth...3 Gel040824_200um_1.csv             200 um         3 0.2427029
## depth...4 Gel040824_200um_2.csv             200 um         1 0.2040657
## depth...5 Gel040824_200um_2.csv             200 um         2 0.2093743
## depth...6 Gel040824_200um_2.csv             200 um         3 0.2081051
## depth...7 Gel040824_200um_3.csv             200 um         1 0.2263164
## depth...8 Gel040824_200um_3.csv             200 um         2 0.2147954
## depth...9 Gel040824_200um_3.csv             200 um         3 0.2215040
##                    Flux Breakpoint Oxygen_Surface Oxygen_breakpoint
## depth...1 -4.644982e-06  -232.9567        265.988           215.260
## depth...2 -4.860979e-06  -229.2693        267.437           214.600
## depth...3 -5.440670e-06  -267.6632        276.688           217.204
## depth...4 -4.574541e-06  -271.7906        265.946           215.835
## depth...5 -4.693543e-06  -268.1471        266.835           215.400
## depth...6 -4.665093e-06  -262.7707        264.796           213.796
## depth...7 -5.073335e-06  -319.0819        282.498           216.533
## depth...8 -4.815069e-06  -236.1913        265.580           213.431
## depth...9 -4.965455e-06  -224.9881        260.208           216.236
Full_results_segmented <- rbind(results_df_1000um_segmented, results_df_200um_segmented)

Full_results_segmented
##                              File Scaffold_Thickness Replicate     slope
## depth...1  Gel040824_1000um_1.csv            1000 um         1 0.5821080
## depth...2  Gel040824_1000um_1.csv            1000 um         2 0.5695080
## depth...3  Gel040824_1000um_1.csv            1000 um         3 0.5511640
## depth...4  Gel040824_1000um_2.csv            1000 um         1 0.5025440
## depth...5  Gel040824_1000um_2.csv            1000 um         2 0.4722500
## depth...6  Gel040824_1000um_2.csv            1000 um         3 0.4822000
## depth...7  Gel040824_1000um_3.csv            1000 um         1 0.3841909
## depth...8  Gel040824_1000um_3.csv            1000 um         2 0.4481669
## depth...9  Gel040824_1000um_3.csv            1000 um         3 0.4936017
## depth...11  Gel040824_200um_1.csv             200 um         1 0.2072080
## depth...21  Gel040824_200um_1.csv             200 um         2 0.2168434
## depth...31  Gel040824_200um_1.csv             200 um         3 0.2427029
## depth...41  Gel040824_200um_2.csv             200 um         1 0.2040657
## depth...51  Gel040824_200um_2.csv             200 um         2 0.2093743
## depth...61  Gel040824_200um_2.csv             200 um         3 0.2081051
## depth...71  Gel040824_200um_3.csv             200 um         1 0.2263164
## depth...81  Gel040824_200um_3.csv             200 um         2 0.2147954
## depth...91  Gel040824_200um_3.csv             200 um         3 0.2215040
##                     Flux Breakpoint Oxygen_Surface Oxygen_breakpoint
## depth...1  -1.304912e-05  -222.7569        350.297           235.901
## depth...2  -1.276666e-05  -214.8453        342.694           230.712
## depth...3  -1.235544e-05  -239.9065        360.905           227.053
## depth...4  -1.126553e-05  -223.8684        327.998           231.279
## depth...5  -1.058643e-05  -217.2484        317.507           225.020
## depth...6  -1.080948e-05  -214.3318        318.232           223.906
## depth...7  -8.612406e-06  -234.9614        312.305           218.707
## depth...8  -1.004656e-05  -231.9015        328.562           218.860
## depth...9  -1.106507e-05  -229.3515        341.538           220.488
## depth...11 -4.644982e-06  -232.9567        265.988           215.260
## depth...21 -4.860979e-06  -229.2693        267.437           214.600
## depth...31 -5.440670e-06  -267.6632        276.688           217.204
## depth...41 -4.574541e-06  -271.7906        265.946           215.835
## depth...51 -4.693543e-06  -268.1471        266.835           215.400
## depth...61 -4.665093e-06  -262.7707        264.796           213.796
## depth...71 -5.073335e-06  -319.0819        282.498           216.533
## depth...81 -4.815069e-06  -236.1913        265.580           213.431
## depth...91 -4.965455e-06  -224.9881        260.208           216.236
mean_flux <- Full_results_segmented %>%
  group_by(Scaffold_Thickness) %>%
  summarize(mean_flux = mean(Flux))

# Plotting the flux values for each thickness
flux_plot <- ggplot(Full_results_segmented, aes(x = Scaffold_Thickness, y = Flux)) +
  geom_boxplot() + geom_point(aes(color = File,shape = Replicate), position = position_jitter(width = 0.05)) + 
  geom_text(data = mean_flux, aes(label = paste("Mean Flux = \n", sprintf("%.2e", mean_flux)), y = mean_flux-.000002), vjust = 0.1, position = position_nudge(x = 0.25)) +
  labs(x = "Thickness (µm)", y = "Flux") +
  ggtitle("Flux Comparison between Thicknesses,\n same volumetric density (4/13/24)")

flux_plot

ggsave("Plots/Polymerized040824_Measured041324_Flux_plot.pdf", flux_plot, device = "pdf")
## Saving 7 x 5 in image
# Running ANOVA
anova_result <- aov(Flux ~ Scaffold_Thickness, data = Full_results_segmented)
summary(anova_result)
##                    Df    Sum Sq   Mean Sq F value   Pr(>F)    
## Scaffold_Thickness  1 1.794e-10 1.794e-10   175.3 4.88e-10 ***
## Residuals          16 1.637e-11 1.020e-12                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# What about for oxygen surface concentration?

mean_Oxygen_Surface <- Full_results_segmented %>%
  group_by(Scaffold_Thickness) %>%
  summarize(mean_Oxygen_Surface = mean(Oxygen_Surface))

# Plotting the surface O2 values for each thickness
ggplot(Full_results_segmented, aes(x = Scaffold_Thickness, y = Oxygen_Surface)) +
  geom_boxplot() + geom_point(aes(color = File), position = position_jitter(width = 0.05)) + 
  geom_text(data = mean_Oxygen_Surface, aes(label = paste("Mean Surface Oxygen = \n", round(mean_Oxygen_Surface, 1)), y = mean_Oxygen_Surface-25), vjust = 0.1, position = position_nudge(x = 0.0)) +
  labs(x = "Thickness (µm)", y = "Oxygen Surface Concentration (µmol/L)") +
  ggtitle("Oxygen Surface Concentration (µmol/L) Comparison between Thicknesses")

# Running ANOVA
anova_result <- aov(Oxygen_Surface ~ Scaffold_Thickness, data = Full_results_segmented)
summary(anova_result)
##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## Scaffold_Thickness  1  18952   18952   119.1 8.03e-09 ***
## Residuals          16   2546     159                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plotting profiles with DBL

files <- list.files(path = "./Polymerized040824_Measured041324", pattern = "\\.csv$", full.names = TRUE)

plot_list <- list()

for (file in files) {
  file_name <- tools::file_path_sans_ext(basename(file))
  
  file_data <- read_csv(paste(file),show_col_types = FALSE) %>%
    mutate(Replicate = as.factor(Replicate)) %>%
    rename("Oxygen (µmol/L)" = `A FireSting Ch1 (µmol/L)`) %>% filter(Condition == "Light" | Condition == "Dark")
  
  # Extract the DBL depth for the current file
  dbl_depth <- Full_results_segmented$Breakpoint[gsub("\\.csv$", "", Full_results_segmented$File) == file_name]

  plot <- ggplot(data = file_data, aes(y = `Oxygen (µmol/L)`, x = `Depth (µm)`,
                                       color = `Condition`, shape = `Replicate`)) + 
      geom_point() + geom_line() +
      coord_flip() + scale_x_reverse() + xlim(0, -1000) + ylim(200, 400) + 
      theme_bw() + ggtitle(file_name) +
      annotate("text", x = -900, y = 395, label = paste("Start:", file_data$Time[1]), vjust = 1, hjust = 1, size = 3) +
      annotate("text", x = -850, y = 395, label = paste("End:", file_data$Time[nrow(file_data)]), vjust = 1, hjust = 1, size = 3) +
      geom_vline(xintercept = dbl_depth, linetype = "dashed") +  # Add DBL depth line
      theme(legend.position = "bottom")

  #print(plot)
  
  plot_list[[file_name]] <- plot
}
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## New names:
## Scale for x is already present. Adding another scale for x, which will replace
## the existing scale.
## • `B No sensor` -> `B No sensor...7`
## • `B No sensor` -> `B No sensor...9`
arranged_plots <- arrangeGrob(grobs = plot_list, ncol = 3)
## Warning: Removed 72 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 72 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 72 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 72 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 72 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 72 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 33 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 33 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 33 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 33 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 33 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 33 rows containing missing values or values outside the scale range
## (`geom_line()`).
gridExtra::grid.arrange(arranged_plots)

ggsave("Plots/Polymerized040824_Measured041324_plots_DBL.pdf", arranged_plots, width = 16, height = 9, units = "in", device = "pdf")